/*!
 * @file main.cpp
 * 
 * Examples of how to use the random number generators.
 * main runs main6 and then main5.
 * main6 shows how to generate various distributions.
 * main5 shows how to chose among various uniform(0,1) generators.
 *
 *
 */
 
#include <cstdlib>
#include <iostream>
#include <iomanip>
#include <string>

#include "mathutil.h"
#include "randomutil.h"
#include "uniform.h"
#include "normal.h"
#include "exponential.h"


// this controls the number of samples printed out in main6.
const int NumDemoMain6 = 10;
// this controls the number of samples printed out in main5.
const int NumDemoMain5 = 10;

// set the number of digits printed out in each number
const int MAIN5_prec=6;
const int MAIN5_width=MAIN5_prec+3;


/* Expected output for main, which runs main6 then main5 using 10 samples.
This is appropriate for importing into matlab to do data visualization and
tests for randomness.


%Running test: main6
%Demo: some random numbers
% U(0,1)  U(0,1)   N(0,1)  N(0,1)   Exp(2) ChiSq(1)
% 0.8643   0.7277  0.2561  0.6094  0.2603  0.0523 ;
% 0.6517   0.7849 -2.6195  0.2186  0.9598  0.0002 ;
% 0.5704   0.5531  0.1316  0.9757  1.1929  0.1478 ;
% 0.0998   0.1129 -0.8087  0.0899  0.1535  0.2599 ;
% 0.6120   0.2960  0.8130 -0.3281  0.0188  3.9259 ;
% 0.7623   0.3977 -0.3762  1.3775  1.1404  0.2986 ;
% 0.2667   0.6739  0.0739  0.7305  0.1779  2.1568 ;
% 0.9776   0.6091 -0.1525 -1.0163  1.1071  2.4748 ;
% 0.6039   0.7382 -0.3435  1.5541  1.1788  0.5500 ;
% 0.1570   0.4420  0.2422  0.4022  0.2757  2.5644 ;

%Done with test: main6

%Running test: main5
% u(0,1)   u(a,b)   N(0,1)   N(mu,sigma),   Exp(lambda) with
a = 1.0000;
b = 4.0000;
mu = 1.0000;
sigma = 0.1000;
lambda = 2.0000;

x = [
 0.505211   2.766053  -1.241855   1.031478   0.165590  ;
 0.428235   3.476248  -0.013056   1.179204   0.067618  ;
 0.652079   1.266737   1.008044   1.143995   3.496783  ;
 0.590803   1.053990  -1.543310   1.091466   0.830395  ;
 0.553513   1.207559  -0.673573   1.024528   0.041253  ;
 0.785994   2.453394  -1.026153   0.983519   0.019392  ;
 0.665093   2.377324   0.589228   0.853047   0.186341  ;
 0.121124   1.982462  -0.549910   1.023855   1.113923  ;
 0.763238   1.833613   0.553358   1.085838   0.137476  ;
 0.086765   1.234408  -0.299389   1.107397   0.301196  ;
];


x_lin = [
 0.613941   1.538193  -0.305221   0.968911   0.362864  ;
 0.699740   3.803533   0.392035   0.971685   0.112633  ;
 0.732603   1.058558   1.144191   1.017488   1.525212  ;
 0.199351   3.251985  -0.853262   1.064979   0.585273  ;
 0.895592   1.108780  -0.082487   0.991642   0.497684  ;
 0.752627   3.027067  -1.402694   1.076453   0.389333  ;
 0.617687   2.253609   0.714094   0.821704   0.277639  ;
 0.302000   2.297700   0.431090   1.039469   0.141527  ;
 0.603310   2.372495   0.885503   0.923203   0.168174  ;
 0.114666   2.351844  -0.615101   1.020700   0.229911  ;
];


x_bound = [
 0.281925   1.845774   0.576057   1.057606   0.165590  ;
 0.191046   1.573139   1.483371   1.148337   0.106007  ;
 0.222326   1.666978  -1.686128   0.831387   0.125724  ;
 0.462464   2.387392   0.405208   1.040521   0.310380  ;
 0.501673   2.505019  -0.996891   0.900311   0.348250  ;
 0.588684   2.766053  -0.621146   0.937885   0.444197  ;
 0.505211   2.515632   0.818528   1.081853   0.351811  ;
 0.126491   1.379473   0.834009   1.083401   0.067618  ;
 0.825416   3.476248  -0.557548   0.944245   0.872675  ;
 0.428235   2.284704   0.269955   1.026996   0.279513  ;
];


*/


/////// first set up some helper functions.
/*!
 * This test driver instantiates the some standard
 * random number generators through the factory methods
 * provided in their respective distributions.
 * It prints the resulting numbers into standard out,
 * which you can pipe to a file or into Matlab.
 */

using namespace std;
using namespace qt::random;

template<typename T1, typename T2, typename T3, typename T4, typename T5>
void testDistributions(T1 & d1, T2 & d2, T3 & d3, T4 & d4, T5 & d5, int n,
                       const string & label) {

    cout << endl;
    cout << label << " = [ " << endl;

    int width = MAIN5_width;
    int prec = MAIN5_prec;

    for (int i=0; i<n; ++i) {
        cout << setw(width) << setprecision(prec) << fixed<< d1() << "  "
             << setw(width) << setprecision(prec) << fixed<< d2() << "  "
             << setw(width) << setprecision(prec) << fixed<< d3() << "  "
             << setw(width) << setprecision(prec) << fixed<< d4() << "  "
             << setw(width) << setprecision(prec) << fixed<< d5() << "  "
             << ";" << endl;
    }

    cout << "];" << endl << endl;

}

/*!
* This routine produces uniform(0,1) samples from a variety of different 
* base generators.  
*/
int main5() {

// this controls the number of samples.
// take n>=500 (or much larger) to do any serious statistical tests
    int n = NumDemoMain5;

    // configure the generator
    unsigned long seed = 38445u;
    double a = 1.0;
    double b = 4.0;
    double lambda = 2;
    double mu = 1;
    double sigma = 0.1;


    // get a RNG
    DefaultGeneratorType & defgen = DefaultGeneratorType::getDefaultGenerator();

    // First we construct default generators for all distributions.
    // The underlying distribution is the mersenne_twister.
    // These all use the same underlying random number generator, defgen.
    // That is, the seed is common to all generators.
    StdUniform::DefaultGenerator ugen = StdUniform::getDefaultGenerator();
    Uniform::DefaultGenerator ugen_ab = Uniform::getDefaultGenerator(a,b);
    StdNormal::DefaultGenerator ngen = StdNormal::getDefaultGenerator();
    Normal::DefaultGenerator ngen_ms = Normal::getDefaultGenerator(mu,sigma);
    Exponential::DefaultGenerator expgen = Exponential::getDefaultGenerator(lambda);


    // Now we create the same generators, except they are bound to
    // the underlying generator invlingen, which is of the type inverse_linear.
    // These generators all use the same
    typedef Generator<inverse_linear> InvLinGenerator;
    InvLinGenerator invlingen(seed);

    StdUniform::GeneratorWrapper<InvLinGenerator>::Generator ugen_invlin
        = StdUniform::getGenerator(invlingen);
    Uniform::GeneratorWrapper<InvLinGenerator>::Generator ugen_ab_invlin =
        Uniform::getGenerator(invlingen,a,b);
    StdNormal::GeneratorWrapper<InvLinGenerator>::Generator ngen_invlin =
        StdNormal::getGenerator(invlingen);
    Normal::GeneratorWrapper<InvLinGenerator>::Generator ngen_ms_invlin =
        Normal::getGenerator(invlingen,mu,sigma);
    Exponential::GeneratorWrapper<InvLinGenerator>::Generator expgen_invlin =
        Exponential::getGenerator(invlingen,lambda);

    // These generators have private seeds that work independently
    // of each other.

    StdUniform::DefaultBoundGenerator ugen_bound
        = StdUniform::getBoundGenerator<DefaultGeneratorType>(seed);
    Uniform::DefaultBoundGenerator ugen_ab_bound
        = Uniform::getBoundGenerator<DefaultGeneratorType>(a,b,seed);
    StdNormal::DefaultBoundGenerator ngen_bound
        = StdNormal::getBoundGenerator<DefaultGeneratorType>(seed);
    Normal::DefaultBoundGenerator ngen_ms_bound
        = Normal::getBoundGenerator<DefaultGeneratorType>(mu,sigma,seed);
    Exponential::DefaultBoundGenerator expgen_bound
        = Exponential::getBoundGenerator<DefaultGeneratorType>(lambda,seed);



    cout << "% u(0,1)   u(a,b)   N(0,1)   N(mu,sigma),   Exp(lambda) with " << endl;
    cout << "a = " << a << ";" << endl;
    cout << "b = " << b <<  ";" << endl;
    cout << "mu = " << mu <<  ";" << endl;
    cout << "sigma = " << sigma << ";" << endl;
    cout << "lambda = " << lambda <<  ";" << endl;


    testDistributions(ugen,ugen_ab,ngen,ngen_ms,expgen,n,string("x"));

    testDistributions(ugen_invlin,ugen_ab_invlin,
                      ngen_invlin,ngen_ms_invlin,expgen_invlin,n,string("x_lin"));

    testDistributions(ugen_bound,ugen_ab_bound,
                      ngen_bound,ngen_ms_bound,expgen_bound,n,string("x_bound"));




}

/*! This class illustrates a simple way
 * to use the library to generate Chi-Squared deviates.
 */

template<class NormalGeneratorType>
class ChiSquaredGenerator {
    NormalGeneratorType & _ngen;
    int _degree;
public:
    /*!
     * @param ngen is passed by reference.  ngen must
     * remain valid while this object is in use.
     */
    ChiSquaredGenerator(NormalGeneratorType & ngen, int degree)
        : _ngen(ngen), _degree(degree) {
        assert(_degree>0);
    }

    /*! adam: Check this implementation.  I may have a constant wrong.  */
    double operator()() {
        double tmp = 0.0;
        double tmp2 = 0.0;
        for(int i = 0; i<_degree; ++i) {
            tmp = _ngen();
            tmp2 += tmp*tmp;
        }
        return tmp2;

    }

};

/*!
* This routine shows different ways to use the random number generators.
*
*/
int main6() {


    // Suppose you don't like the default generators.
    // Now we mix and match uniform generators with
    // methods for generating other distributions.

    unsigned long seed = 264383u;

    typedef  Generator<minstd_rand> MinGeneratorType;
    typedef  Generator<mersenne_twister2> MersenneGeneratorType;

    // instances of some base random number generators
    MinGeneratorType minGenerator(seed);
    MersenneGeneratorType mersenneGenerator(seed);

    // here's one way to roll your own normals using a different base generator.
    PolarStdNormalGenerator
    <MinGeneratorType > polargen(minGenerator);


    // Here we use the static version of StdNormal::invCdf.
    // We could use any function of type double (*)(double).
    // We use the
    // generate exponential deviates using the InvCdfGenerator class.
    InvCdfGenerator
    <std::pointer_to_unary_function<double,double>, MersenneGeneratorType>
    invgen(std::ptr_fun(StdNormal::invCdf), mersenneGenerator);

    // Now roll your own exponentials.

    // Here we use the nonstatic versions of Exponential::invCdf, lambda = 2.0
    // It's not pretty, but it's very flexible.
    // If you want pretty, use the factory methods.
    double lambda = 2.0;
    Exponential e(lambda);
    std::mem_fun1_t<double,qt::random::Exponential,double> memfun =
        std::mem_fun(&Exponential::invCdf);

    InvCdfGenerator
    <std::binder1st
    <std::mem_fun1_t<double,qt::random::Exponential,double> >,
    MinGeneratorType >
    invgen2( std::bind1st(memfun,&e), minGenerator);

    // Now Chi-Squared
    ChiSquaredGenerator<PolarStdNormalGenerator <MinGeneratorType > >
    chisqgen(polargen,1);



    cout << "%Demo: some random numbers" << endl;
    cout << "% U(0,1)  U(0,1)   N(0,1)  N(0,1)   Exp(2) ChiSq(1) " << endl;
    for (int i=0; i<NumDemoMain6; ++i) {
        cout <<"%" << setw(7) << setprecision(4) << fixed<< minGenerator() << "  "
             << setw(7) << mersenneGenerator() << " "
             << setw(7) << polargen() << " "
             << setw(7) << invgen() << " "
             << setw(7) << invgen2() << " "
             << setw(7)  << chisqgen()  << " ;" << endl;
    }
}

int main(int argc, char *argv[])
{

    cout << endl << "%Running test: main6 " << endl;
    main6();
    cout << endl << "%Done with test: main6 " << endl;
    cout << endl << "%Running test: main5 " << endl;
    main5();
    cout << endl << "%Done with test: main5 " << endl;

    // comment out this line if you want to pipe the output into a MATLAB .m script.
    // keep it if you use the command line and want it to stay open.
    // system("PAUSE");

    return EXIT_SUCCESS;
}
